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Abstract 

The fundamental higher-order Landau plasma modes are known to be generally heavily damped. 



We show that these modes for the ion component in a weakly ionized plasma can be substantially 
modified by ion-neutral collisions and a dc electric field driving ion flow so that some of them can 
become unstable. This instability is expected to naturally occur in presheaths of gas discharges at 
sufficiently small pressures and thus affect sheaths and discharge structures. 
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I. INTRODUCTION 



One of the fundamental kinetic phenomena in plasma physics is the higher-order Landau 
modes (l|. They are the heavily damped solutions of the dispersion relation describing 
the electrostatic modes of a one-component collisionless Maxwellian plasma [2|. (A one- 
component plasma is the approximation where only one plasma component oscillates.) The 
dispersion relation is: 
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is the plasma dispersion function, U is any contour in the complex £ plane passing from 
£ = —00 to +00 below the singular point £ = x, u is the complex wave frequency, k is the 



real wave number, v t = \Jk-gT / n is the thermal velocity, T is the temperature, \i is the 
particle mass, k-Q is the Boltzmann constant, Ad = i>t/wp is the Debye length, and u p is the 
plasma frequency. The transcendental equation (TJQ) with respect to uj yields the Langmuir 
mode and an infinite number of the heavily damped higher modes |l|. 

Because the higher modes are a fundamental phenomenon, it is not surprising that they 
received considerable attention in the literature despite their strong damping. The first, but 
implicit, mention dates back to Landau himself [2], who used terms "all the poles" and "that 
of the poles" in relation to the solutions of Eq. ([1]). An explicit statement on the existence 
of the higher modes was made 14 years later by Jackson, who demonstrated their presence 
aua.ytica.lv & Nu,„e r ica. r esu. tS were pub.ished in the 1960s Q fl Q. Receat.y, these 
modes have been studied in relativistic plasmas p, |7J. Experimentally, the higher modes 
were observed in the 1970s (in the spatially damped case) [8|, [9|. 

Very similar, but nonlinear, kinetic electron waves have been recently observed in simu- 
lations to persist without any apparent decay over many plasma periods after being excited 
by an artificial external driver |10l.lll|. They were called KEEN waves, standing for "kinetic 



electrostatic electron nonlinear" waves 
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llj . One remarkable similarity between them 



and the higher modes is that KEEN waves were obtained even for frequencies well below 
the plasma frequency and that their typical phase velocity was somewhat above the electron 
thermal velocity, which is in full accord with the properties of the higher modes. Another 
similarity can be seen in the quasineutral character of KEEN waves (evidenced by Fig. 9 of 



Ref. jll|). as the higher modes are quasineutral at small wave numbers [i.e., the neglect of 
the vacuum term in Eq. (JT|) does not affect the higher solutions u at small wave numbers, 
as can be easily verified]. The heavy Landau damping, characteristic to the higher modes, 
was eliminated for KEEN waves by a population of trapped particles created by the exter- 
nal drive. Very recently, analogous nonlinear kinetic waves have been studied for the ion 



component in the presence of the electron response using a similar external driver [12] . 

In this paper, we show that a finite number of the higher modes can become unstable 
under a quite natural set of circumstances, so that kinetic waves can emerge without exci- 
tation by an artificial external driver. Namely, we show that the higher modes for the ion 
component in a weakly ionized plasma can be substantially modified by ion-neutral collisions 
and ion flow driven by a dc electric field so that some of these modes can become unstable. 



Such flows are an essential feature of low-pressure gas discharges [13M17J , as they naturally 
arise to maintain the balance of absorption of ions and absorption of electrons on the walls 
and electrodes of the discharge chamber. 

So far, there have been numerous investigations of streaming instabilities triggered due to 
relative flows of various plasma components in col 



isionless plasmas, with perhaps the most 
known example being the Buneman instability [18], but our study is principally different in 
two important aspects explained below. 

The first aspect is that we include a dc electric field and collisions with neutrals, and do 
this self- consistently. The self-consistency here means two things. First, we find the steady- 
state velocity distribution from the model itself, i.e., from the balance of ion acceleration in 
the field and ion- neutral collisions (instead of assuming a model distribution, e.g., a shifted 
Maxwellian distribution). Second, collisions and the dc field not only define the steady state 
but are also fully accounted for in the analysis of perturbations (this is essential to obtain the 



instability, as we show in Sec. lIHC]) . Of course, such an approach generally requires extremely 



cumbersome velocity calculations, but we avoid this difficulty by considering the common 
case where the dominant mechanism of collisions is charge transfer, in which the ion and 
neutra , Sim p, y exchange identic and thus velocity BQ. Our further approbation 



is to assume the collision frequency to be velocity independent, which al 
solution not only for the steady state but also for the ion susceptibility 



ows an elegant 
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22J. As, in 



reality, it is not the collision frequency but the cross section that is characterized by a 



weak (logarithmic) velocity dependence 20 



23] (in the regime where charge transfer is the 



dominant mechanism of collisions), we separately show that the instability remains in the 
constant mean free path case (in Sec. IIIIDj) . 

The second aspect is that we consider only one oscillating plasma species, i.e., ions. The 
electron density is assumed to be fixed, which physically implies that the temperature of 
the electrons is high enough so that they do not "feel" electric fields, similar to how ion 
Langmuir waves are derived in classical textbooks |24|. This assumption allows us to render 
a simple physical picture, while we address the role of the electron temperature in Sees. IIIIEI 
and IIV Al (with additional details given in Appendix [C]) and provide explicit conditions for 
neglecting the electron effects [Eqs. (|35|) . ( 136|) . (1C4j) . and ( 1C5j) ]. As our model involves 
only one oscillating plasma component, it is remarkable indeed that our analysis reveals an 
instability. This instability is clearly associated with a novel mechanism. 

We emphasize that the ion-kinetic instability described in this paper is not a variation 



of any known collisionless instability (e.g., the bump-on-tail instability 25|, |26|]), since a 
collisionless one-component plasma with our steady-state velocity distribution is always 
stable (Sec. IIIICj) . We explain the instability mechanism in Sec. IIIIC1 

This ion-kinetic instability should affect a large class of gas discharges, as it is expected 

^29] under quite common conditions (Sec. IIVBI) . 



to naturally occur in presheaths [14|, 



15. 
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In even broader context, our study shows that the often disregarded higher-order Landau 
modes can in fact play a crucial role in the presence of an electric field. 

II. METHODS 

A. Basic equations 

Let us consider a weakly ionized plasma in a dc electric field E driving ion flow. For 
electrons we assume a Boltzmann distribution with a sufficiently large temperature so that 
we can consider their number density Uq to be homogeneous and fixed (note that Boltzmann 
electron distributions in the presence of field-driven ion flow are common in discharges 



14J . ll5l.l27H29|; we address the role of the electron temperature in Sees. IHIEI and llVAj) . We 



assume Eq to be homogeneous and use the kinetic equation for ions with the Bhatnagar- 



Gross-Krook (BGK) ion-neutral collision term 21] and Poisson's equation: 



?f + v .^ + ±/ r Eo-^V — 

dt dr m \ dr J dv 

= -vf + u$> M jf(v')dv' 1 (3) 

-A(j> = ^(Jfdv-n Y (4) 

where / is the ion distribution function, <fr is the electric potential describing the time-space 
varying field (i.e., the field apart from E ), 

$ m(v)= /0 _.,. 2v v? ex p(-^r) (5) 



(2tt<) 3 / 2 ^ V 2w 



is the normalized Maxwellian velocity distribution of neutrals, v is the velocity-independent 
ion-neutral collision frequency, v tn = yk^T^Jm is the thermal velocity of neutrals, T n is the 
temperature of neutrals, e is the elementary charge (ions are assumed to be singly ionized), 
m is the ion mass, and eo is the permittivity of free space. The BGK term exactly describes 
charge transfer collisions under the assumption of a velocity-independent collision frequency 



21] 



B. Steady state 

The homogeneous steady-state solution f = f is found from Eqs. ([3]) and (j4j) by setting 
= 0, df/dt = 0, df/dr = 0. This gives [2l|, bflj: 



A(v> = (s5f r exp ( _? - jr if ! ) * <6) 

where 

eE 
v f = (7) 

mv 

and the subscript "f" stands for "flow," as the flow velocity (1/no) / v /o<^ v can be shown 



to be equal to Vf. It is helpful to note that Eq. (EJ) here is simply Eq. (3) of Ref. (211 ]. 
but rewritten using another integration variable in order to show that /o is an integral 
superposition of shifted Maxwellian distributions with exponential weights. At E = 0, the 
velocity distribution (jSJ) is Maxwellian, with the thermal velocity being equal to that of 



J 0,z tn 




FIG. 1: Steady-state solution. Shown is the longitudinal velocity distribution fo >z = 
f_ J_ fodv x dv y , where the z axis is in the direction of the field. Different lines correspond 
to different values of the parameter u = Vf/vt n . 

neutrals. The solution ([6]) is shown in Fig. [TJ As Fig. [1] shows, at large fields the solution (pQ) 
cannot be approximated by a shifted Maxwellian distribution, as the solution fl6]) becomes 
highly asymmetric with respect to the position of its maximum. This can also be illustrated 
by considering the limit of cold neutrals, v tn — ¥ 0, where Eq. fl6]) becomes 



/o(v) = — exp ( — - j S(v x )S{v y ), v z >0, 

/o(v) = 0, v z < 0, 



(8) 



where the z axis is in the direction of Eq. 



C. Dispersion relation 



The dispersion relation is obtained by linearizing Eqs. ([3]) and (jl]) wit 



i respect to 4> and 



/ — /o and can be written using the ion susceptibility derived in Refs. [21|, |3l| . A derivation 
of this dispersion relation by solving the initial value problem is given in Ref. 22) . The 
dispersion relation is: 
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where Wpi = ^/noe^J{e^m} is the ion plasma frequency, w is the complex wave frequency, and 
k is the real wave number; the solutions u of Eq. (Q provide contributions oc exp(— iut+ik-r) 
to the asymptotic expression for the solution <fi = 0(r, t) of the initial value problem at large 



t |22| 



The result <^ is different from what one obtains by simply substituting our steady- 
state distribution (jBJ) to the classical expression for the dielectric function of a collisionless 



plasma 24]. The difference is due to our accounting for the perturbation term (eE /m) ■ 
d(f — fo)/dv and the perturbation of the right-hand side of Eq. (J3J). It is this difference 
that results in the instability, as shown in Sec. IIIICI 

D. Analysis 



First, we numerically analyze the dispersion relation (|9]) (in Sec. IIII A|) . The corresponding 
dimensionless variables are the flow parameter 



u=— , (10) 

Utn 



the collision parameter 



v 



C= -' (U) 

the dimensionless frequency uj/uj p i and the dimensionless wave number kA, where 



A = — . (12) 

w pi 

Note that in the above definition of A we use the thermal velocity of neutrals, so at Eq = 
the length A becomes the ion Debye length. For finite Eq, the ion thermal velocity and 
the ion Debye length are not defined because of the non-Maxwellian form of the velocity 



distribution. Also note that the length A is not the effective screening length in the presence 



of the flow, as evident from Ref. 32J. The dimensionless form of the dispersion relation (Q 
is given in Appendix |A] 

Second, we provide an analytical proof of the instability existence using Eq. §§§ (in 
Sec. IIIIBp and explain the instability mechanism (in Sec. IIIICp . 

Third, we analyze whether the instability remains in the constant mean free path case 



(in Sec. IIIIDj) . For this purpose, we replace the right-hand side of Eq. (J3J) by [33] : 



St[/(r,v)] = J J^i [$M(v)/(r,vO - $ M (v')/(r, v)] dv', (13) 

where £ is the collision length. This operator exactly describes charge transfer collisions 
under the assumption of a velocity-independent cross section. 

Finally, we numerically study the effect of the electron response (in Sec. IIII El) . To clarify, 



there are two effects related to a finite electron temperature: (i) the electron response to ion 
oscillations, and (ii) a finite inhomogeneity length of the Boltzmann electron distribution in 
the field E . Here we focus on effect (i) by adding the Boltzmann response term l/(/cA e ) 2 
to the left-hand side of Eq. ( I9a| . where A e = [eo&B^e/ (^oe 2 )] 1 ^ 2 is the electron Debye length 
and T e is the electron temperature. Effect (ii) is discussed in Sec. IIV Al 

III. RESULTS 

Let us briefly summarize our results before describing them in detail: 

1. The instability occurs when u > 8 and C, < 0.3 (Sec. IIII Al) . However, in the limit 
( — > 0, u = const [physically corresponding to a collisionless one-component plasma 
with our steady-state velocity distribution fl6])] the instability growth rate tends to 
zero (Sec. [THAI Sec- HTTC]) . 

2. It is downstream waves that become unstable at the above instability boundary 
(Sec. IIII A]) . To clarify, by "downstream waves" we mean that the phase velocity vector 
Re(o;)k/fc 2 is in the direction of E . 

3. When the above instability conditions are met, the instability occurs in a finite range 
of wave numbers. Its lower and upper ends for downstream waves can be estimated 
as k ~ u/vf and k ~ (z//t> f )min {w 2 , C~ 4//3 }, respectively [Sec. IIII Cj . 

8 



4. The number of unstable modes depends on u and ( and can in principle be made 
arbitrarily high (Sec. IIIIBj) . 



5. The frequency of the most unstable mode for downstream waves is a; ~ (3.35 + 
O.G4i)^eE k/m, provided that the above instability conditions on u, ( and k are 
satisfied by a considerable margin (Sees. IHI Bl and IIIICI) . 

6. The physics of the resulting growing waves is essentially kinetic (Sec. IIIICI) . 

7. The instability remains in the constant mean free path case (Sec. IIIIDI) . 

8. The electron response for typical values of T e does not shift the above instability 
thresholds to unrealistic values, as for T c /T n = 200 (corresponding to k-oT e = 5 eV, 
T n = 300 K) the first higher mode remains unstable at, for instance, u = 14 and 
C = 0.05 (SecHHED- 

Let us now describe these findings in detail. 

A. Numerical results 

This section provides the results of the numerical analysis of the dispersion relation ([9]) . 

No-flow case: At u = and ( — > 0, Eq. (jHJ) is equivalent to the Landau dispersion relation 
(JT]). The solutions in this case are shown in the left column of Fig. [2j Note that in the limit 
k — y 0, the higher modes are acoustic, i.e., u oc k, with the proportionality coefficients being 
complex numbers with comparable real and imaginary parts |l|. 

A finite ( merely results in that Im(o;) for any given mode (including the ion Langmuir 
one) tends to a constant at k — > 0, as shown in the right column of Fig. |2j This constant for 
all higher modes is the same and equal to — v. For the ion Langmuir mode, this constant 
differs by a factor of two and is equal to —u/2, for v < 2cu pi . 

Effect of field: As a detailed discussion of our numerical results in light of the instability 
mechanism is given in Sec. IIIICI in this section we only provide a brief description of what 
happens in the presence of the field. 

Let us first consider the dispersion curves for downstream waves. The left column of 
Fig. |3] shows that at u = 2 and ( = 0.1, the ion Langmuir mode is no longer the least 
damped mode at large wave numbers, as the first two higher modes have smaller decay 
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FIG. 2: Modes in the absence of flow. Shown are the solutions of Eq. ([9]) for u = 0. The left 
column represents the collisionless case [£ — > 0; in this case the dispersion relation is reduced to 
Eq. ([ID], and the right column illustrates the effect of collisions for £ = 0.1. The ion Langmuir 
mode and the first two higher modes are denoted by "L," "1," and "2," respectively. The dashed 
lines, shown to guide the eye, correspond to Re(w)//c = v tn . 

rates. (We determine which mode is the ion Langmuir mode by looking at the behavior at 
k —?• 0: the ion Langmuir mode is the one that has a finite real part of the frequency at 
k — y 0.) The right column of Fig. [3] illustrates that at u = 10 and £ = 0.1, the first higher 
mode is unstable in a range of wave numbers, while the ion Langmuir mode remains stable 
at these parameter values. It is not only the first higher mode that can become unstable, 
as we found numerically a large number of unstable modes by increasing u and decreasing 
( (in accordance with the analytical results of Sees. IIII B[) . 

For upstream waves (for which the phase velocity vector is in the direction opposite to the 
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FIG. 3: Modes in the presence of flow. The left and right columns correspond to u = 2 and u = 10, 
respectively, both are for £ = 0.1 and downstream waves. The notation of the modes is the same as 
in Fig. The right column illustrates the instability of the first higher mode. The upper and lower 
dash lines, shown to guide the eye, correspond to Re(u)/k = V{ and Ke(ui)/k = v tn , respectively. 

flow) we did not find any instability. The general case of arbitrary angle of propagation can 
be mathematically reduced to the case of propagation along or against the flow, as Eq. (jUJ) 
contains Vf only in the combination k ■ Vf. 

We calculated the instability region in the (u,() space (see Fig. @J. One can see that 
the instability region is bound within u > 8 and Q < 0.3. Here, two important comments 
need to be made. First, we found numerically that in the limit ( — > and u = const within 
the instability region, the dimensionless growth rate tends to zero, Im(u/uj p [) — > + . Thus, 
a finite collision parameter is essential for the instability. Second, it is always downstream 
waves that become unstable at the boundary of the instability region. 
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FIG. 4: Stability diagram. The instability region is bound within u > 8 and £ < 0.3. Note that in 
the limit £ — > 0, u = const within the instability region, the growth rate of the instability tends to 
zero, so that a finite collision parameter is essential for the instability. 

B. Analytical proof of the instability existence 



Let us analytically prove the existence of the instability and show that the number of 
unstable modes can be made arbitrarily high by varying parameters u and £• We consider 
Eq. (Q in the limit v tn — > 0, w P i — )■ oo at finite Vf, is, u, k and assume that k is in the 
direction of the flow. Then, the first term (unity) in Eq. fl9a|) is negligible so that the 
dispersion relation takes the form B = 0. Furthermore, in Eq. (19dl) the second term of the 
expression inside the square brackets is negligible as well. In the resulting dispersion relation, 
let us consider the limit of large k. This allows us to neglect the unity in the denominator 
of Eq. (|9cj) as well as the unity in the first term of the right-hand side of Eq. (l9dl) and yields 



u> 



C/k 



VflS, 



where the numerical factor C is given by 



exp iC£ - -i? ) dt, = 



(14a) 



(14b) 



Equation ( 114bj) has an infinite number of solutions C, and they all have positive real and 
imaginary parts. That is, we get an infinite number of unstable modes. They correspond 
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to downstream waves, as Re(w) > for all solutions. The most unstable solution, i.e., the 
one with the largest Im(C), is C ~ 3.35 + 0.64i. The presence of an infinite number of 
unstable modes in the above limit means that, by approaching this limit, one can achieve an 
arbitrarily large number of unstable modes. This is consistent with our numerical analysis 
of the dispersion relation OH]), as we numerically found a large number of unstable modes at 
sufficiently large u and sufficiently small ( (as noted in Sec. lIH A|) . Also, we have numerically 
verified that the solutions described by Eq. ( 114)1 indeed satisfy Eq. (Q in the limit considered 
above. We also provide a derivation of Eq. ( 1T4")) starting from our initial equations flSJ)-©], 
see Sec. IHTU 

C. Instability mechanism 

First of all, let us emphasize that our instability is not a variation of any known col- 



lisionless instability (e.g., the bump-on-tail instability [25|, |26|), as we show below that a 
collisionless one-component plasma with our steady-state velocity distribution (JSJ) is always 
stable. We consider Eq. (Q in the limit of infinitely small Eo and u, but keep their ratio 
(which determines Vf and hence the form of the steady-state velocity distribution) finite. In 
this limit Eq. simplifies to 

1 r ie W {i[u/{kv tn )]i-e/2}di _ 
^ (kXyJ l + ^(k-v f )/(^ tn ) ~ l ' 

This equation does not have unstable solutions, as we verified numerically. Note that this is 
in accordance with the above stated fact that in the limit ( — > 0, u = const the instability 
growth rate tends to zero. 

To explain the instability mechanism, let us start with the physics of the higher modes 
in the case of u = and £ — > and then see how they can become unstable in the presence 
of the field and collisions. 

For u = and ( — > 0, the higher modes represent quasineutral oscillations at small wave 
numbers, as mentioned in the Introduction and explained in more detail below. Let us first 
clarify that the term "quasineutral" here means that the neglect of the A0 term in Poisson's 
equation does not affect the higher modes, or, to put it differently, that contributions from 
various velocity domains to the velocity integral determining the susceptibility compensate 
for each other. To demonstrate this, we first note that the A<p term in Poisson's equation 

13 



contributes to Eq. ([T]) by means of the first term of Eq. ([I]) (unity), as seen from the 
derivation of Ref. |2j. This vacuum term, as seen directly from Eq. (JTJ, indeed does not 
affect the higher-order solutions u in the limit k — Y 0, because of the l/A; 2 -factor in the 
second term of Eq. ([I]) and the acoustic character of the higher modes {to oc k) at k — > 0. 
Appendix iBlexplicitly shows how contributions from various velocity domains to the velocity 
integral determining the susceptibility compensate for each other. 

These quasineutral waves are heavily Landau damped, as their phase velocity is of the 

order of the thermal velocity. One way to eliminate the damping is, as noted in the In- 

;roduction, to apply an external driver in order to create a population of trapped particles 
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12] . Another way, which is what is considered in this paper, is to induce field-driven 
flow. 

Let us now explain how growing quasineutral ion waves become possible in the presence 
of field-driven ion flow, starting from our initial equations ([3])-(j5]). 

We consider the case where the flow velocity is much larger than the thermal velocity 
of neutrals and focus on the kinetics of ions with velocities in the range v tn < ^ < ff. 
Mathematically, this is equivalent to considering the limit of cold neutrals and simplifying 
the resulting steady-state distribution ([8]) by taking its low velocity part, 

/o(v) = —8(v x )6(v y ), {v t »H(> Vtn) > 0, 

Vf 

/o(v) = 0, v z < 0. (16) 

By using this expression and considering downstream waves with f — fo = / a exp(— itot+ikz), 
(ft = (p^exp^icot+ikz) (where the subscript "a" stands for "amplitude"), we get the following 
linearized kinetic equation: 

eE df a ikecp^riQ 

- %ujj & + ikv z f & H d(v) 

m ov z m Vf 



■vh + v&(v) j UWdV. (17) 



Note that for this consideration to be valid, an additional term 
[ikecj) au n /{mv f)] exp(—v z /vf)S(v x )5(v y ) (for v z > 0) should be negligible, as the latter 
appears when Eq. (JS} is used instead of Eq. ( JIB]) . We will analyze this condition, as well as 
other restrictive conditions imposed below, later in this section. 

Let us consider the ballistic case, where the right-hand side of Eq. ( 1T7J1 is negligible. This 
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implies 



and 



v 



/adv 



< 



fce|0 a |ra o 



(18) 



f <C max {| u; |, kv z } , (19) 

which should be satisfied for the characteristic velocity v z of the solution / a ; Eq. (TIB"]) follows 
from comparison of the terms containing the delta-function [i.e., the fourth and sixths terms 
in Eq. (fTTl)]. and Eq. ( Fl9|) follows from comparison of the terms not containing the delta- 
function [i.e., the first, second, and fifth terms in Eq. HIT)) : we did not insert the third term 
of Eq. (ITT]) into Eqs. ( TT8|) and ( TT9|) . because this term is the only remaining term, so that its 
magnitude is determined by the terms already included in Eqs. ( ITB"]) and (119))]. The solution 
of Eq. (ITTj) with zero right-hand side is 



/a 



z/c0 a n 



exp 



m 



lUJV, 



eE 
x6(v x )S(v y ), 



ikv\ 
~Y 

v z > 



By substituting this solution to Poisson's equation k 2 
sion relation 



U = 0, v z < 0. (20) 

= (e/e ) / fa, dv, we get the disper- 



1- 



e E kv 



exp 



f Jo 



zcutv — 



ifcf^ 



^2 = 0, 



(21) 



_eE V 2 

where the main contribution to the integral should be from the aforementioned range v tn "C 
v z ^C Vf [see Eq. ( 1161) ] for the consideration to be valid. Note that the vacuum term in 
Eq. (I2ip is the first term (unity), as it comes from the fc 2 a -term in Poisson's equation. 

Let us consider the case where neglecting the vacuum term in Eq. (|21j) does not affect 
some of its solutions u (the condition for this is obtained below). Neglecting the vacuum 
term, we get 

(22) 



u = C 



m 



where C is given by Eq. (114b)) . Thus, we obtained exactly Eq. ( 1X4]) [Eqs. (122)) and (I14ap are 
equivalent] and hence its infinite set of unstable solutions. The above consideration explains 
the instability mechanism, as we identified the elementary terms and processes essential for 
the instability and how they work together. 

Let us now show how the restrictive conditions imposed above determine the instability 
region in the (u, £, k) space. Namely, in the following, we show that for the most unstable 
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mode the above restrictive conditions can be summarized as: 

l«^«min{u 2 ,C 4/3 }. (23) 

One can immediately see from condition (|23|) that it can only be satisfied when 

u » 1, ( < 1, (24) 

which explains why the instability occurs at large u and small ( (see Fig. H]). As condition 
( 123]) is for neglecting the effects that are not essential for the instability mechanism, the 
instability range can be estimated as 

l<^<min{V,C 4/3 }. (25) 

To derive Eq. (123]) . let us start with the condition for neglecting the vacuum term. We 
rewrite Eq. ( ]2"T1) as 

3/2 



'-'(^"fcff)' 



x / exp ( j - j - 1— ] d^ = 0, (26) 

Jo \ v^l^M 2 / 

and note that for the most unstable mode 

Re(w) ~ Im(w) ~ ^eE k/m. (27) 

It follows that the vacuum term does not affect the most unstable mode when 



v J yv,k) » L (28) 

Let us now consider the condition that the characteristic velocity v z providing the main 
contribution to the integral in Eq. ( 121]) satisfies v tn <?; z < Vf. By using Eqs. ( 127]) we obtain 
that this velocity is v z ~ *JeEo/(mk), which gives 

kVr 9 . . 

K — < u 2 . (29) 

Concerning the remaining conditions — Eqs. ( 1T8]) and ( 1T9]) and the condition of smallness 
of the term mentioned just after Eq. ( I17j) - - they are automatically satisfied for the most 
unstable mode when Eqs. ( 128]) and (129]) hold, as can be easily shown. 
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By combining Eqs. (|28|) - (|29|) we get Eq. ( )23|) . Note that another derivation of Eq. ( 123|) 
[from Eq. (Q] is provided in Appendix O 

Let us now discuss, in light of the above findings, our numerical results shown in the right 
column of Fig. El To do so, we first rewrite Eqs. ( 12~2"]) and ( 12"5"]) in the dimensionless units 
used to plot Fig. [3J 



ijj 



— = C^/u(k\ (30) 

^<k\<mm\(u,-^\, (31) 

where, to remind, C ~ 3.35 + 0.64i for the most unstable mode. One immediately sees that 
the real part of the frequency shown in the right column of Fig. [3] is in excellent agreement 
with Eq. (!30l) (inside the instability range), as the difference is within 15%. Remarkably, 
the square root dependence of Re(o;) on k is easily noticeable on the graph. The imaginary 
part, however, does not show good agreement, which is not surprising because the parameter 
values corresponding to the right column of Fig. [3]are close to the instability boundary shown 
in Fig. H] [Eq. (|22l) is derived assuming that all effects that are not essential to the instability 
mechanism are negligible]. To verify this interpretation of the discrepancy, we run a test 
comparing Im(o;) given by Eq. (130]) with Im(o;) given by Eq. (Q for parameters u, ( and kX 
very well satisfying inequalities ( l3Tj) and found excellent agreement. Finally, comparison of 
Eq. ( 13~T]) with the right column of Fig. [3] shows reasonable agreement, as Eq. ( 131]) for the 
parameters of the right column of Fig. [3] becomes 0.01 < kX < 0.2. 

Let us also briefly discuss the wave number corresponding to the maximum growth rate as 
well as the phase velocity of the most unstable mode. Concerning the former, Eq. ( 122]) shows 
that the growth rate increases with k, so the wave number corresponding to the maximum 
growth rate can be estimated as the upper end of the instability wave number range (125]) . 
As regards the phase velocity of the most unstable mode, we combine Eqs. (122]) and ( 1251) 
and thus obtain that this phase velocity varies with k in the range 

max{^ n ,< 2 / 3 }<^<^. (32) 

Finally, let us emphasize that the unstable solutions (12"2~]) are quasineutral. This already 
follows from the derivation itself, as these solutions are obtained for the case where they 
are not affected by the neglect of the vacuum term (see above). Let us, however, explicitly 
illustrate how high- and low-velocity particle contributions to the susceptibility compensate 
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for each other. The susceptibility is the second term in Eq. (|2Tj) . and we consider it as the 
sum of two parts, one being the contribution from v z < Ke(cj)/k and one from v z > Ke(u)/k. 
By using the most unstable solution uj pa (3.35 + O.GA^y/eEok/m, we obtain that the first 
part normalized by (tu P i/ u) 2 [u / (v{k)] 3j/2 is ~ 0.1 + 0.04i and that the second part normalized 
as above is ~ —0.1 — 0.04z, so that they add up to zero. 

D. Constant mean free path case 

It is already clear from Sec. IIIICI that the instability mechanism is generic in the sense 
that it is not sensitive to a particular velocity dependence for charge transfer collisions. Let 
us, however, explicitly demonstrate that the instability remains in the constant mean free 



path case. In the limit of cold neutrals the operator ffT3|) simplifies to 32] 

St[f} = - V -l + 6 -^lJf( r ,v')v'dV. (33) 

Then the steady-state distribution is 

/o(v) = — - exp z — 8(v x )S(v y ), v z > 0, 

/o(v) = 0, v z <0, (34) 

where Vf t i = \ fvfodv\/no = ^2eE £/(7im) is the flow velocity in the constant mean free 
path case. The subsequent steps are exactly the same as in Sec. IIIICI and we come to the 
conclusion that the relation ( 1221) exactly applies to the constant mean free path case as well 
(in the limit considered). 

E. Role of the electron response 

Let us numerically show that the electron response for typical values of T e does not shift 
the instability thresholds to unrealistic values. (How the electron response is included is 
explained in Sec. Ill PI ) We choose a typical value T e /T n = 200 corresponding to /c B T e = 5 eV 
and T n = 300 K. Performing calculations for u = 10 and ( = 0.1 (i.e., for the same values 
of u and ( as in the right column of Fig. |3]), we find that the instability disappears (see 
Fig. |5]), but the instability does not disappear at, for instance, u = 14 and ( = 0.05 (see 
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0.02 



0.01- 



-0.01 



-0.02 




FIG. 5: Effect of the electron response, (a) No response case: The unstable mode for u = 10, 
£ = 0.1 (i.e., for the same u and £ as in the right column of Fig. [3]). (b) The same mode for 
T c /T n = 200, with u and £ as above, (c) The same mode for u = 14, £ = 0.05 and T e /T n as above. 

Fig. EJ). These values of u and ( are quite realistic, as will be seen in light of the discussion 
of Sec. HVBl 

ig. It should also 



24j requires a finite 



We did not find any other unstable mode in cases (a), (b) and (c) of ! 
be noted that the classical expression for the ion-acoustic instability 
electron-to-ion mass ratio for the instability to occur, while in our approach the electron 
mass is effectively zero. 



IV. DISCUSSION 



A. Applicability limits 



The applicability of our model is limited by the following factors: 

(i) electron response to ion oscillations (addressed in Sec. IIIIEI) . 

(ii) inhomogeneity of the Boltzmann electron density profile in the field Eo (discussed 
below), 

(iii) assumption of a velocity-independent collision frequency for charge transfer collisions 
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(addressed in Sec. IIIIPp . and 



(iv) assumption that charge transfer is the dominant mechanism of ion scattering on neu- 
trals (discussed below). 

Concerning factor (iv), for argon at room temperature, charge transfer indeed dominates 
when ion velocities exceed the thermal velocity of neutrals by a factor of ~ 3 or larger 34]. 
Thus, for argon at room temperature, factor (iv) should not substantially influence our insta- 
bility, as the latter occurs at u > 8. As regards factor (ii), the corresponding inhomogeneity 
length is L e = k&T e /(eE ). Our model can be applied when this distance is larger than both 
the ion-neutral collision length which is Vf/u (for Vf > v tQ ) and the wavelength 2ir/k. This 
imposes the following requirement: 

T c (,«!/] 

— > max <^ u , — - \ . (35) 

T n { v tn k J 

Let us make estimates to find when the instability occurs in the case of a finite elec- 
tron temperature. To do so, we need to find when the instability mechanism explained in 
Sec. IHICl is not substantially influenced by other effects, including the electron effects, in a 
range of wave numbers. This means that 

• condition (TSU) is satisfied and 



• in at least a part of the wave number range (|23l) both condition (I35jl is met and 
the effect of the Boltzmann response term l/(/cA e ) 2 on the most unstable mode is 
negligible. 

As shown in Appendix [Cj conditions (a) and (b) above can be written as 

w»l, C<1, 7^>« 2 , (36) 

[note that the last inequality in Eq. ( l3T)|) coincides with one of the two conditions imposed 
by Eq. ( )35l) ]. Here two comments need to be made. First, the fact that the last inequality of 
Eq. (1361) is better satisfied for the stable curve (b) than for the unstable curve (c) of Fig. [5] 
does not mean a contradiction between Eq. ( 136]) and Fig. [5j as the first two inequalities of 
Eq. ( )36|) are better satisfied for the curve (c) than for the curve (b). Second, estimates ( 136|) 
are for neglecting the effects that are not essential to the instability mechanism, so the 
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instability threshold on the electron temperature can be estimated from Eq. (|36|) as 

Y ~ « 2 (37) 

(for a given u ^> 1 and a given (<1). Thus, as the instability thresholds on u and ( in 
the absence of the electron effects are about 8 and 0.3, respectively (as given by the exact 
calculation presented in Fig. HJ), we can summarize the instability conditions as 

8<u<M, C<0.3; (38) 

V -^n 



Eq. (1381) is not very accurate (as seen from Fig. [5]) because the condition on T e in Eq. (138]) 
is merely an estimate and not an exact result. 

B. Presheaths 

Based on the above, we expect the instability to naturally occur in presheaths at suffi- 
ciently small, but still quite common, pressures, as explained in the following. Let us see 
how the instability conditions ( |38l) can be met in presheaths. 

We first analyze the condition £ < 0.3. By replacing v by ffcn n , where a is the ion- 
neutral cross section and n n is the neutral number density, we can rewrite the condition 
C < 0.3 as 

where P is the gas pressure (here we took the condition u > 8 into account, i.e., we sub- 
stituted Vf ~ 8f tI1 , because at larger Vf the restriction on P is stronger). By substituting 
T n = 300 K and a = 6.5 x 10~ 15 cm 2 (this value of a is derived in Ref. [32] for argon from 
the data of Ref. [35]), we get 



p ^ 2P Wl5^ (40) 

The obtained condition can be easily satisfied in gas discharges, as there have been many 
experiments with pressures below 2 Pa and plasma densities about or greater than 10 14 m~ 3 



36H40J. (Note that here n denotes the local density in the presheath and not the density in 



the bulk of the discharge, but 



on the presheath structure 



or a rough estimate one can use the latter; see measurements 



[liU. 



) 
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Let us now discuss the condition 8 < u < \JT C /T n in light of the Bohm criterion [2j 
First of all, we note that the Bohm criterion applies when the ion-neutral collision length is 
larger than the electron Debye length [28]. This condition can be written as 



P< 



1 n e 2 k B T£ 



a 



eoT e 



(41) 



This is a weaker condition than Eq. (139]) because they differ by the factor 30 yT n /T e , which 
is typically larger than unity, so the Bohm criterion applies when Eq. (139]) is satisfied. 
According to the Bohm criterion, the ion flow velocity reaches at least the Bohm speed 
\lk B T c jm at the sheath-presheath edge [27] . Then, the condition 8 < u < ^JT C /T n is met 
in a certain space region within the presheath if 




> 

1 n 



(42) 



For T n = 300 K, this means ksT c > 1.6 eV, which is often satisfied. 



The instability may have significant imp 



tant to plasma physics and technology 20 



ications, as presheaths and sheaths are impor- 



29J. First, the instability may result in a flow 



turbulence or various dynamic structures and thus lead to the appearance of strong electric 
fields varying on the ion time scale. An alternative is the formation of a static structure that 
suppresses the instability. In an extreme scenario, the instability may significantly affect the 
whole discharge or even switch it off. The above effects cannot be modeled using the hydro- 
dynamic approach since the latter ignores the higher modes. (Inaccuracy of hydrodynamic 
modeling of plasma boundary layers is illustrated in, e.g., Ref. |4l|.) 



C. Presheath instability experiments 



Instabilities in presheaths have been observed in the presence of two ion species 14j, |42 |. 
The ion-kinetic instability described in this paper, in contrast, can occur even for a single 
ion species plasma. Interestingly, measurements of Ref. 42] indicate that the presheath 
instability did not disappear when only one ion species was present (although in this case 
the instability amplitude was reduced; see Fi g. 1 of Ref. |42j). To explain this, the presence 
of unidentified impurity ions was suggested [42j. We suggest that these results might also 
be explained by the instability of the higher modes. 
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Let us explicitly suggest experimental conditions to observe the instability described in 
this paper. Concerning the choice of gas, argon is an excellent candidate, as in this gas at 
room temperature, charge transfer is indeed the dominant mechanism of ion scattering; when 
ion velocities exceed the thermal velocity of neutrals by a factor of ~ 3 or larger |34|. In 
addition, an argon plasma well represents a single ion species plasma, as the impurity degree 
due to ions other than Ar + (e.g., Ar ++ and Ar^~) is usually quite small 43). Concerning the 
gas pressure, it should satisfy Eq. fj4*0l) . An example of an experiment well satisfying this 
condition is Ref. |44j. [Note, however, that at extremely small pressures, the instability may 
not be detectable, as its growth rate tends to zero at ( — > (Sees. UTTAI and ILTlCp .] It may be 
also necessary to take measures to increase the electron temperature, as our estimate of the 
instability threshold on the electron temperature yields k-QT e ~ 1.6 eV [Eq. (l4"2"j) ] so that the 
actual threshold may differ by a factor of a few. These measures can be: (i) decreasing the 
gas pressur e |36| . (ii) decreasing the rf peak-to-peak voltage [36|, and (iii) using a Maxwell 



demon 
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D. Dusty plasmas 



Another implication of the instabi 



affect the interparticle interaction 
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ity described in this paper is that this instability 



46[ and ion drag force [21] in dusty plasmas 47H50]. 



may 



In particular, our results imply that the expression for the shielding of a dust particle in 
the presence of ion flow given by Eq. (6) of Ref. 32] is only valid when the ratio of the 
"field-induced Debye length" (defined in Ref. 32() to the collision length is larger than a 
certain threshold, which is supposedly close to 0.3, i.e., to that in the BGK case (otherwise, 
the linear response formalism does not apply because of the instability of the steady state). 
The results of Ref. [46[ are unaffected because, in that work, the subthermal flow regime 
was considered. The resulting change in the interaction between dust particles may affect 
their self-organization and dynamics 51n62]. 



V. CONCLUSION 

We found a remarkable type of instability, which can be triggered in a weakly ionized 
plasma in the presence of ion flow driven by a dc electric field. We showed that the dc field 
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and ion-neutral collisions can substantially modify the physics of the ion higher-order Landau 
modes so that some of them can become unstable. The instability is of broad relevance to 
gas discharge physics, as dc fields are common in gas discharges. In particular, the instability 
is expected to naturally occur in presheaths under quite common conditions [Eqs. (T4"Uj) and 
(Pi2")) ] and thus affect sheaths and discharge structures. In even broader context, our study 
shows that the often disregarded higher-order Landau modes can in fact play a crucial role 
in the presence of an electric field. 

Acknowledgments 

The authors thank the Editorial Board Member for a thorough review and helpful com- 
ments. R.K. acknowledges the support of the Science Foundation for Physics within the 
University of Sydney and the hospitality of the Max-Planck-Institut fur extraterrestrische 
Physik, Germany. A.V.I, and G.E.M. received funding from the European Research Council 
under the European Unions Seventh Framework Programme (FP7/2007-2013) / ERC Grant 
agreement 267499. The work was partially supported by the Australian Research Council. 

Appendix A: Dimensionless form of Eq. (|5j) 

The dispersion relation fl9]) can be written in our dimensionless units for a numerical 
analysis as 

i+ <^( i -£f"*<-*>«r 

x p«p(-gK 

J 1 + iuk, cos 6 

where 9 is the angle between k and Eo- 

Appendix B: Quasineutrality of the higher modes 

This Appendix illustrates, for u = and ( — > 0, the quasineutral character of the higher 
modes at small wave numbers. Namely, we explicitly show how contributions from various 
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velocity domains to the velocity integral determining the susceptibility compensate for each 
other. The susceptibility, i.e., the second term in Eq. (JTJ, can be written as J2j: 

e " f dfazjvz) dv z 

me k 2 J dv z v z — u/k 1 

u 

where the wave number is in the direction of the z axis, fo,z( v z) — 
(27r) _1 / 2 (no/ftn) exp[— Vg/(2v^ a )] is the velocity distribution integrated over v x and v y , 
and U is any contour in the complex v z plane passing from v z = — oo to +00 below the 
singular point v z = u/k. This integral can be written as the sum of the integral on the 
real axis and the contribution from the pole of the integrand [for Im(u) < 0] . Dividing the 
integral on the real axis into two parts, one over v z smaller than the phase velocity R,e(u)/k 
and one over the remaining interval, we can write 

X = j^(Ji + J2 + h), (B2) 



where 



, < [ Re(uj)/k dhM) dv z 

n J-oo dv z v z - u/k 

T v? Q [°° df 0>z (v z ) dv z 

^0 jRe(co)/k dv z V z - U/k 



2Tciv? n df , z (v z ) 
J?. — 



n dv z 



(B5) 

v z =ui/k 



the integration in the expressions for J\ and J2 is performed on the real axis and J3 is the 
contribution from the pole. By substituting here the first higher mode solution u/(kv tn ) w 
3.60 - 1.73i (at k -> 0), we get J x w -0.032 - 0.06z, J 2 ~ 5 x 10" 5 - 3 x 10 _4 z, J 3 ~ 
0.032 + 0.06i, so that it is explicitly seen that they add up to zero and that, as a side note, 
J2 is almost unimportant here. For the second mode, we have u/(kv tn ) ~ 4.47 — 2.87i, 
Jx w -0.012 - 0.035i, J 2 « 5 x 10" 7 - 6 x 10" 6 z, J 3 w 0.012 + 0.035i. 

Appendix C: Estimates 

This appendix provides estimates for the instability region and thus derives Eqs. ( 123]) 
and (136]) . We first consider the "pure" case where all terms that are not essential for the 
instability are neglected. This allows us to obtain estimates for the magnitudes of the 

25 



essential terms. We then compare them with the neglected terms to obtain conditions under 
which the neglected terms are indeed negligible. 

The "pure" case can be considered by omitting all terms in Eq. (^ that were neglected in 
Sec. IIIIBI These terms are: (i) the unity in round brackets in Eq. (I9dj) . (ii) the second term 
of the expression inside the square brackets in Eq. (19dj) . (iii) the unity in the denominator 
in Eq. (!9c|) . and (iv) the first term (unity) in Eq. (I9a| . Therefore, the main contribution to 
the integral in Eq. f )9c|) is from 

To derive this, we expressed to via k using Eq. (TTIj) ; we assume Re(C) ~ Im(C) ~ 1, as 
we consider the most unstable mode; we also assume that k is in the direction of Eo- The 
above estimate of rj yields 

/ n 3/2 

' B '~fe) ' (C2) 

This result is obtained by replacing 77 and drj by the estimate ( 1C1|) and substituting unity 
for the exponent. Analogously, for the second term in Eq. fl9a|) we get 

B w 2 

,2 



^pi 



.4 



uk 3 v 3 




^ minn, \l — }. (C3) 



'f 

Let us now make a comparison with the magnitudes of the non-essential terms (i)-(iv). 
The term (i) is negligible when \u\ ^> v. This is equivalent to k 3> u/vf. The term (ii) 
can be omitted when k ^C V{u/v 2 n . The term (iii) does not play any role when it is smaller 
than the other term in the denominator of Eq. (|9cl) with rj replaced by estimate ( 1C1I) . This 
gives k 3> v/i)f, which coincides with the condition for the neglect of the term (i). Finally, 
the term (iv) can be neglected when it is smaller than the right-hand side of Eq. ( 1C3|) . The 
right-hand side of Eq. (1C3|) can be simplified using the condition k ^> u/i)f for the neglect 
of the term (i). The result is that the term (iv) is negligible when k <ti (z//ff)(w P i/z/) 4 ' 3 . By 
combining the above conditions, we get Eq. (1231) . 

Analogously, we obtain that the electron response term l/(/cA e ) 2 , added to the left-hand 
side of Eq. fl9al) . is unimportant when: 

vfum 2 

k>> Jw- (C4) 

Let us now see when conditions (1C4I) and (|35|) are met in at least a part of the wave number 
range (1231 assuming that condition ff24|) is satisfied. As conditions ( 1C4|) and fl35|) impose 

2(3 



lower (and not upper) limits on the wave number, we only need to find when conditions ( 1C4I) 
and (133|) are met at the upper end of the wave number range (123]) . Concerning Eq. ( IC4j) , it 
is met at the upper end of the wave number range (j23p when 

^»max{n,n 2 C 2 / 3 }. (C5) 

■*■ n 

As regards Eq. (1331) . it is met at the upper end of the wave number range ([23]) when 
T e /T n ^> w 2 . This inequality is a stronger condition than Eq. ( 1C5I) . Thus, we arrive at 
Eq. fl33J. 
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